Step 5 Phylogenetic Trait Graph

About

Workspace

#Install packages if needed

remotes::install_github("eliocamp/ggnewscale@dev")

remotes::install_github("YuLab-SMU/ggtree")

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("Biostrings")
BiocManager::install("ggtree")

#For document rendering
tinytex::install_tinytex()

# general
library(tidyverse)
library(readxl)
library(plyr)
library(dplyr)
library(here)
"%notin%" = Negate("%in%")
here::here()
library(tinytex)

# graphics
library(ggplot2)
library(ggnewscale)
library(viridis)

# phylogenetic & other tools
library(ape)
library(Biostrings)
library(ggtree)
library(rotl)  #for phylogenetic analyses, get all the species? from Hinchliff et al. 2015 PNAS
library(stringr)
library(tidytree)

Load Prey & Traits Data

NOTES:

• Using combination of my trial code and Matt Savoca and company’s code to build a phylo tree with associated trait and %FO data.

• A lot of troubleshooting later and I found that the species list needs to go on the far left side of the data, all other data/factors need to be added to the RHS of the species list

# data contains sp list, class, order & family + traits.

# Probable life stage data
my_prey_prob = read.csv(here("./data/output_data/prey_trait_select.csv"), header = TRUE) %>%
    dplyr::select(prey_sp, prey_class:prey_family, life_stage:maxM, -c(diet_likely,
        life_note)) %>%
    filter(prey_sp %notin% c("Lampanyctus mexicanus", "Janthina exigua"))
# Need to remove these two species as they cause issues with ROTL due to naming

# 306 species and 13 variables of interest excluding prey_sp str(my_prey_prob)

Most of the variables except the diet metrics need to be factors.

sapply(my_prey_prob, class)

##            prey_sp         prey_class         prey_order        prey_family 
##        "character"        "character"        "character"        "character" 
##         life_stage       vert_habitat       horz_habitat       diel_migrant 
##        "character"        "character"        "character"          "numeric" 
##   diel_migrant_cat         refuge_cat     season_migrant         season_cat 
##        "character"        "character"          "numeric"        "character" 
## gregarious_primary              maxFO               maxN               maxM 
##        "character"          "numeric"          "numeric"          "numeric"

my_prey_prob[, 2:13] <- lapply(my_prey_prob[, 2:13], factor)

my_prey_prob$maxFO[my_prey_prob$maxFO == 0] <- NA
my_prey_prob$maxN[my_prey_prob$maxN == 0] <- NA
my_prey_prob$maxM[my_prey_prob$maxM == 0] <- NA

Tree Build

Building the basic tree

# Tree build 1
breaks <- c(seq(1, nrow(my_prey_prob), 50), nrow(my_prey_prob) + 1)  # why are we doing this?
# looking up each of the species in dataset to a phylogeny, doing it by sets of
# 50 if there's an NA (species that didn't match a value in known phylogeny),
# breaks the whole chunk of 50 species

for (i in 1:(length(breaks) - 1)) {
    taxa <- as.character(my_prey_prob$prey_sp[breaks[i]:(breaks[i + 1] - 1)])
    taxa <- taxa[taxa != "" & !is.na(taxa)]

    resolved_namest <- tnrs_match_names(taxa)  # I think this is where all the extra species fall out
    resolved_namest <- resolved_namest[!is.na(resolved_namest$unique_name), ]  # ignore an NA
    if (i == 1) {
        resolved_namess <- resolved_namest
    } else {
        resolved_namess <- rbind(resolved_namess, resolved_namest)
    }
}
resolved_names <- resolved_namess
resolved_names <- resolved_names[resolved_names$flags != "INCERTAE_SEDIS_INHERITED",
    ]

# original tree based on simple taxon datset
my_tree_prob <- tol_induced_subtree(ott_ids = resolved_names$ott_id, label_format = "name")

my_tree_prob$tip.label <- gsub("_", " ", my_tree_prob$tip.label)  # removes underscore between genus and species names
my_tree_prob$tip.label <- str_extract(my_tree_prob$tip.label, "[A-Z][a-z]+ [a-z]+")

my_tree_prob <- compute.brlen(my_tree_prob, method = "Grafen", power = 1/2)  #add branch lengths to my tree using the Grafen (1989) method
my_tree_prob <- ladderize(my_tree_prob, right = TRUE)

View(my_tree_prob)

Save Tree

write.tree(my_tree_prob, here("./data/output_data/albacore_diet_tree_prob"))

Load Saved Tree

my_tree_prob <- read.tree(here("./data/output_data/albacore_diet_tree_prob"))
#306 species

Sort out species name issues for graphing

NOTE: Here we edit the tree tip labels - species names - for my_tree_prob and check data dimensions line up with original df containing trait info my_prey_prob.

# Edit tip labels
my_tree_prob$tip.label <- as.factor(sub("_", " ", my_tree_prob$tip.label))
# Check that the number of rows matches the number of tip labels
nrow(my_prey_prob) == length(my_tree_prob$tip.label)  #TRUE

## [1] TRUE

# Can check this manually also nrow(my_prey_prob) # 306
# length(my_tree_prob$tip.label) # 306

NOTE: We have in the past encountered a data frame dimension issue when plotting our data values onto the tree. There are also several species names that do not line up. This issue is addressed in the chunk below.

# my_prey_prob$prey_sp <- gsub(' ', '_', my_prey_prob$prey_sp) #for ease of
# plotting take away ' '
tree_names_prob = data.frame(sort(my_tree_prob$tip.label))  #sort the species names from the phylo tree
# There are 306 sort brings tree_names_prob from 301 to 300 because of an NA in
# the tip labels
prey_names_prob = data.frame(sort(my_prey_prob$prey_sp))  #sort the species names from the original list

nrow(tree_names_prob)

## [1] 306

nrow(prey_names_prob)  #Needs to be 306 -- 306

## [1] 306

# Sometimes one species didn't make it to the tree

NOTE: Here we identify which species names do not match from our original species and trait df, compared to that which parsed through phylogenetic info from ROTL project. These were:

• Original data my_prey_prob, and checked via prey_names_prob: Diplodus sargus, Leuroglossus stilbius, Liocranchia reinhardtii, Mullus barbatus barbatus, Notoscopelus kroyeri, Robustosergia robusta.

• Tree data my_tree_prob, and checked via tree_names_prob: Bathylagus stilbius, Liocranchia reinhardti, Mullus barbatus, NA, Notoscopelus elongatus, Sergia robusta.

tree_names_prob == prey_names_prob

## here we check the names that are in the prey data but not on the tip labels
## of the tree
names_not_in_tree_prob = prey_names_prob %>%
    filter(sort.my_prey_prob.prey_sp. %notin% tree_names_prob$sort.my_tree_prob.tip.label.)
print(names_not_in_tree_prob$sort.my_prey_prob.prey_sp.)

## [1] "Diplodus sargus"          "Leuroglossus stilbius"   
## [3] "Liocranchia reinhardtii"  "Mullus barbatus barbatus"
## [5] "Notoscopelus kroyeri"     "Robustosergia robusta"

## check names that are in the tree tip labels but not in the prey data
names_not_in_prey_prob = tree_names_prob %>%
    filter(sort.my_tree_prob.tip.label. %notin% prey_names_prob$sort.my_prey_prob.prey_sp.)

print(names_not_in_prey_prob$sort.my_tree_prob.tip.label.)

## [1] Bathylagus stilbius    Liocranchia reinhardti Mullus barbatus       
## [4] NA                     Notoscopelus elongatus Sergia robusta        
## 306 Levels: Abralia redfieldi Abraliopsis affinis ... Walvisteuthis rancureli

NOTE: These taxa need to be re-aligned via detailed code below.

## here we decide to match the names in the prey data to the tip labels, so we create a dataset with the names in the prey data that are in the tree (I see this looks kind of confusing with the double negative, the names not not in the tree). These are the names that do not need to be fixed, which we will bind the fixed names to afterwards so we do not get doubles of the same names (with different spelling).
my_prey_keep_prob = my_prey_prob %>%
  filter(prey_sp %notin% names_not_in_tree_prob$sort.my_prey_prob.prey_sp.) 
#300 names to keep, and 6 need to be corrected

## these are the names we are going to fix to match the tip labels of the tree, given they are synonyms/misspellings. These are the same names in 'names_not_in_tree_adult', except they now have the columns of data which will allow use to merge it to 'my_prey_keep_adult' after we fix the names.
my_prey_fix_prob = my_prey_prob %>%
  filter(prey_sp %in% names_not_in_tree_prob$sort.my_prey_prob.prey_sp.)
#6 names to fix

## this step isn't necessary, but it does make it easier in the names are in alphabetical order
my_prey_fix_prob = my_prey_fix_prob[order(my_prey_fix_prob$prey_sp),]
# the names have to be characters for the name reassignment to work, it will give an error if we don't do this (they are factors before this)
my_prey_fix_prob$prey_sp <- as.character(my_prey_fix_prob$prey_sp)
## we make these into x2 and y2 because it makes it much easier to write the next section
x2 = my_prey_fix_prob
y2 = as.vector(names_not_in_prey_prob$sort.my_tree_prob.tip.label.)

unique(x2$prey_sp)

## [1] "Diplodus sargus"          "Leuroglossus stilbius"   
## [3] "Liocranchia reinhardtii"  "Mullus barbatus barbatus"
## [5] "Notoscopelus kroyeri"     "Robustosergia robusta"

y2

## [1] "Bathylagus stilbius"    "Liocranchia reinhardti" "Mullus barbatus"       
## [4] "NA"                     "Notoscopelus elongatus" "Sergia robusta"

#here we look at the names in both of the lists and see if there are misspellings/synonyms in the names that caused the discrepancy between the name lists. As you can see this is true for all of the species in these lists expect one. 'Diplodus sargus' does not have a match in the tree, since the only other tip labels left in the tree is an NA after the rest of the reassignments.
# if you want to see how this works you can run individual pieces of this before running the whole thing. For example x2[4,1] is "Leuroglossus stilbius", and y2[1] is "Bathylagus stilbius". Then after running this code we match the name in the prey data (x2) to the names in the tree (y2), so "Leuroglossus stilbius" becomes "Bathylagus stilbius".

x2[1,1] = y2[4]; #--> Diplodus sargus = NA
x2[2,1] = y2[1]; #--> Leuroglossus stilbius = Bathylagus stilbius
x2[3,1] = y2[2]; #--> Liocranchia reinhardtii = Liocranchia reinhardti
x2[4,1] = y2[3]; #--> Mullus barbatus barbatus = Mullus barbatus
x2[5,1] = y2[5]; #--> Notoscopelus kroyeri = Notoscopelus elongatus
x2[6,1] = y2[6]; #--> Robustosergia robusta = Sergia robusta


# we need to do the opposite for one of the names here "Diplodus sargus", which is in the prey data but not in the tip labels. So, we are replacing the NA in the tip labels
tip_labels_prob <- as.vector(my_tree_prob$tip.label)
tip_labels_prob[28] ##this is the NA we want to replace

## [1] "NA"

tip_labels_prob[28] = "Diplodus sargus"

## replacing the tip labels with the additional name
my_tree_prob$tip.label <- tip_labels_prob #Still have 306 tip labels

# here we bind the fixed names to the names that didnt need to be fixed
my_prey_prob = rbind(my_prey_keep_prob, x2)

#here we get rid of any of the names in the data that didn't have a match in the tree tip labels, which was only 'Diplodus sargus'
  #--> This is where one species drops out!
# Do not run this bit
#my_prey_prob = my_prey_prob %>% 
  #--> This is where one species drops out!
#  filter(prey_sp %in% my_tree_prob$tip.label)

nrow(my_prey_prob);length(my_tree_prob$tip.label) ## one of those tip labels is an NA

## [1] 306

## [1] 306

#the prey names need to be the rownames for the heatmap to work on the tree
#rownames(my_prey_prob) = NULL #this did not solve the row names thing
rownames(my_prey_prob) = my_prey_prob$prey_sp

#tip labels need to be characters to pipe onto the tree
my_tree_prob$tip.label <- as.character(my_tree_prob$tip.label)
str(my_tree_prob)

## List of 5
##  $ edge       : int [1:578, 1:2] 307 308 309 310 311 312 313 314 315 316 ...
##  $ edge.length: num [1:578] 0.00164 0.24959 0.0066 0.01789 0.03479 ...
##  $ Nnode      : int 273
##  $ node.label : chr [1:273] "mrcaott42ott150" "mrcaott42ott49" "mrcaott42ott658" "Clupeocephala" ...
##  $ tip.label  : chr [1:306] "Sebastes wilsoni" "Sebastes zacentrus" "Sebastes proriger" "Sebastes brevispinis" ...
##  - attr(*, "class")= chr "phylo"
##  - attr(*, "order")= chr "cladewise"

Paper figures

Basic phylogenies

Basic prey species phylogeny graph

# Phylotree without labels
prob_basic = ggtree(my_tree_prob, layout = "circular")
prob_basic

# Phylotree with species labels
prob_basic_lab = ggtree(my_tree_prob, layout = "circular") + geom_tiplab(size = 2)
prob_basic_lab

Class-based phylogeny graph

#Use basic + prey_sp labels + prey_class

# creating a dataset that splits the species data by class
prey_class_prob_info <- split(x = my_prey_prob$prey_sp, f = my_prey_prob$prey_class)
unique(my_prey_prob$prey_class)

## [1] Actinopterygii Cephalopoda    Gastropoda     Hexanauplia    Hydrozoa      
## [6] Malacostraca   Thaliacea     
## 7 Levels: Actinopterygii Cephalopoda Gastropoda Hexanauplia ... Thaliacea

# splitting the tree species by class
my_tree_prob_class <- groupOTU(my_tree_prob, prey_class_prob_info)
as.character(sort(unique(my_prey_prob$prey_class)))

## [1] "Actinopterygii" "Cephalopoda"    "Gastropoda"     "Hexanauplia"   
## [5] "Hydrozoa"       "Malacostraca"   "Thaliacea"

#creating a palette for the class split
#pal_prob_class = c("grey40", '#0047ab', '#751308', '#4B0082', '#F05E23', '#013220', '#B80F0A', "#66C2A5", "#3288BD") ## the grey40 is needed for a branch between branches, but doesn't assign to a class


pal_prob_class = c('#0047ab', #Actinopterygii
                    #'#751308', #Branchiopoda #older teal colour "#66C2A5" #Not in prob data just was used for adult data
                    '#4B0082', #Cephalopoda
                    "#AD15E2", #Gastropoda #old orange '#F05E23'
                    '#B80F0A', #Hexanauplia
                    "#773405", #Hydrozoa #'#013220'
                    "#E86103", #Malacostraca '
                    "#3288BD",  #Thaliacea
                   "grey40" ## the grey40 is needed for a branch between branches, but doesn't assign to a class
                    
                    )

## plotting the tree without tip labels

prob_class <- ggtree(my_tree_prob_class, aes(color = group), layout = 'circular') +
   theme(legend.position = c(0.5,0.50),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_colour_manual('Class', aesthetics = c('colour'), values = pal_prob_class,
                      breaks = c("Actinopterygii", "Cephalopoda","Gastropoda", "Hexanauplia", "Hydrozoa", "Malacostraca","Thaliacea"),
                      labels = c("Actinopterygii", "Cephalopoda","Gastropoda", "Hexanauplia", "Hydrozoa", "Malacostraca","Thaliacea"))

prob_class

# plotting the tree with tip labels
prob_class_lab <- ggtree(my_tree_prob_class, aes(color = group), layout = 'circular') +
   theme(legend.position = c(0.5,0.50),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_colour_manual('Class', aesthetics = c('colour', 'fill'), values = pal_prob_class,
                      breaks = c("Actinopterygii", "Cephalopoda","Gastropoda", "Hexanauplia", "Hydrozoa", "Malacostraca","Thaliacea"),
                      labels = c("Actinopterygii", "Cephalopoda","Gastropoda", "Hexanauplia", "Hydrozoa", "Malacostraca","Thaliacea")) +
  geom_tiplab(size = 2)

prob_class_lab

Diet use data

Frequency Occurrence

#Overlay maxFO (frequency of occurrence data) on a basic prey tree coloured by Class
#Use prey_class, maxFO

prob_maxfo <- gheatmap(prob_class, my_prey_prob[,"maxFO", drop = FALSE], 
                       offset = 0, width = 0.05, colnames = FALSE) +
  theme(#legend.position = c(0.5,0.50),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_fill_viridis_c(name = "Percent frequency (%FO)", na.value = 'grey')

prob_maxfo

Abundance

#Overlay maxN (abundance data) on a basic prey tree coloured by Class
#Use prey_class, maxN

my_prey_prob$maxN <- as.numeric(my_prey_prob$maxN)
prob_maxn <- gheatmap(prob_class, my_prey_prob[,"maxN", drop = FALSE], 
                      offset = 0, width = 0.05, colnames = FALSE) +
  theme(#legend.position = c(0.5,0.50),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_fill_viridis_c(name = "Percent abundance (%N)", na.value = 'grey')
prob_maxn

Biomass

#Overlay maxM (biomass data) on a basic prey tree coloured by Class
#Use prey_class, maxM

prob_maxm <- gheatmap(prob_class, my_prey_prob[,"maxM", drop = FALSE], 
                      offset = 0, width = 0.05, colnames = FALSE) +
  theme(#legend.position = c(0.5,0.50),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_fill_viridis_c(name = "Percent biomass (%KG)", na.value = 'grey')
prob_maxm

Trait data

Habitat Use Traits

For habitat use traits: we are going to plot multiple rings of data around the phylogenetic trees using vert_habitat, horz_habitat, diel_migrant (change to yes/no), and season_migrant (yes/no)

#creating this plot one iteration at a time. Adding one layer of the heatmap, then allowing for another legend on the plot, and then another layer of the heatmap

habitat_traits_prob1 <- gheatmap(prob_basic, my_prey_prob[, 'vert_habitat',drop=FALSE], 
                                 offset= 0, width=0.05, colnames = FALSE) +
  scale_fill_viridis_d(name = "Vertical Habitat (VH)", 
                       direction = -1, 
                       breaks = c("benthic", "demersal", "epipelagic", "mesopelagic", "bathypelagic"),
                       limits = c("benthic", "demersal", "epipelagic", "mesopelagic", "bathypelagic"), 
                       na.translate = TRUE, 
                       option = 'D', 
                       na.value = 'grey')+
  theme(legend.position = c(1,0.80),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_y_continuous(expand = c(0,4))+ #this opens up the gap for the Vertical habitat label around the tree
  annotate('text', x = 1.03, y = -7, label = 'VH', angle = -85, size = 4)

habitat_traits_prob1.5 <- habitat_traits_prob1 + new_scale_fill()

habitat_traits_prob2 <- gheatmap(habitat_traits_prob1.5, my_prey_prob[ ,'horz_habitat',drop=FALSE], 
                                 offset=0.05, width=0.05, colnames = F) +
  scale_fill_viridis_d(name = "Horizontal Habitat (HH)", 
                       option = "D", 
                       direction = -1,
                       breaks = c("reef-associated", "coastal", "continental shelf","continental slope", "oceanic"),
                       limits = c("reef-associated", "coastal", "continental shelf","continental slope", "oceanic"),
                       na.translate = TRUE, 
                       na.value = 'grey')+
  theme(legend.position = c(1,0.683),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_y_continuous(expand = c(0,4))+
annotate('text', x = 1.08, y = -7.1, label = 'HH', angle = -85, size = 4)

habitat_traits_prob2.5 <- habitat_traits_prob2 + new_scale_fill()

habitat_traits_prob3 <- gheatmap(habitat_traits_prob2.5, my_prey_prob[ , 'diel_migrant',drop=FALSE], offset=0.10, width=0.05, colnames = F) +
  scale_fill_manual(name = "Diel Migrant (DM)",
                    breaks = c("0", "1"),
                    limits = c("0", "1"),
                    labels = c("no","yes"),
                    values = c("0"="#FDE725FF", "1"="#39568CFF"), 
                    na.value = 'grey')+
  theme(legend.position = c(1,0.5985),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_y_continuous(expand = c(0,4))+
  annotate('text', x = 1.13, y = -7.5, label = 'DM', angle = -85, size = 4)

habitat_traits_prob3.5 <- habitat_traits_prob3 + new_scale_fill()

habitat_traits_prob_final <- gheatmap(habitat_traits_prob3.5, 
                                      my_prey_prob[ ,'season_migrant',drop=FALSE],
                                     offset=0.15, width=0.05, colnames = F) +
  theme(legend.position = c(1.05,0.50),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_fill_manual(name = "Seasonal Migrant (SM)",
                    breaks = c("0", "1"),
                    limits = c("0", "1"),
                    labels = c("no", "yes"),
                    values = c("0"="#FDE725FF", "1"="#39568CFF"), 
                    na.value = 'grey') +
  scale_y_continuous(expand = c(0,4))+
  annotate('text', x = 1.18, y = -7.5, label = 'SM', angle = -85, size = 4)

  
#Check graph
habitat_traits_prob_final

Habitat use and aggregation trait

Adding aggregation to habitat use figure

#creating this plot one iteration at a time. Adding one layer of the heatmap, then allowing for another legend on the plot, and then another layer of the heatmap

habitat_traits_prob1 <- gheatmap(prob_basic, my_prey_prob[, 'vert_habitat',drop=FALSE], 
                                 offset= 0, width=0.05, colnames = FALSE) +
  scale_fill_viridis_d(name = "Vertical Habitat (VH)", 
                       direction = -1, 
                       breaks = c("benthic", "demersal", "epipelagic", "mesopelagic", "bathypelagic"),
                       limits = c("benthic", "demersal", "epipelagic", "mesopelagic", "bathypelagic"), 
                       na.translate = TRUE, 
                       option = 'D', 
                       na.value = 'grey')+
  theme(legend.position = c(1,0.80),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_y_continuous(expand = c(0,4))+ #this opens up the gap for the Vertical habitat label around the tree
  annotate('text', x = 1.03, y = -7, label = 'VH', angle = -85, size = 4)

habitat_traits_prob1.5 <- habitat_traits_prob1 + new_scale_fill()

habitat_traits_prob2 <- gheatmap(habitat_traits_prob1.5, my_prey_prob[ ,'horz_habitat',drop=FALSE], 
                                 offset=0.05, width=0.05, colnames = F) +
  scale_fill_viridis_d(name = "Horizontal Habitat (HH)", 
                       option = "D", 
                       direction = -1,
                       breaks = c("reef-associated", "coastal", "continental shelf","continental slope", "oceanic"),
                       limits = c("reef-associated", "coastal", "continental shelf","continental slope", "oceanic"),
                       na.translate = TRUE, 
                       na.value = 'grey')+
  theme(legend.position = c(1,0.683),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_y_continuous(expand = c(0,4))+
annotate('text', x = 1.08, y = -7.1, label = 'HH', angle = -85, size = 4)

habitat_traits_prob2.5 <- habitat_traits_prob2 + new_scale_fill()

habitat_traits_prob3 <- gheatmap(habitat_traits_prob2.5, my_prey_prob[ , 'diel_migrant',drop=FALSE], offset=0.10, width=0.05, colnames = F) +
  scale_fill_manual(name = "Diel Migrant (DM)",
                    breaks = c("0", "1"),
                    limits = c("0", "1"),
                    labels = c("no","yes"),
                    values = c("0"="#FDE725FF", "1"="#39568CFF"), 
                    na.value = 'grey')+
  theme(legend.position = c(1,0.5985),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_y_continuous(expand = c(0,4))+
  annotate('text', x = 1.13, y = -7.5, label = 'DM', angle = -85, size = 4)

habitat_traits_prob3.5 <- habitat_traits_prob3 + new_scale_fill()

habitat_traits_prob4 <- gheatmap(habitat_traits_prob3.5, 
                                      my_prey_prob[ ,'season_migrant',drop=FALSE],
                                     offset=0.15, width=0.05, colnames = F) +
  theme(legend.position = c(1.05,0.50),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_fill_manual(name = "Seasonal Migrant (SM)",
                    breaks = c("0", "1"),
                    limits = c("0", "1"),
                    labels = c("no", "yes"),
                    values = c("0"="#FDE725FF", "1"="#39568CFF"), 
                    na.value = 'grey') +
  scale_y_continuous(expand = c(0,4))+
  annotate('text', x = 1.18, y = -7.5, label = 'SM', angle = -85, size = 4)

habitat_traits_prob4.5 <- habitat_traits_prob4 + new_scale_fill()

habitat_traits_prob_agg_final <- gheatmap(habitat_traits_prob4.5, 
                                      my_prey_prob[ ,"gregarious_primary",drop=FALSE],
                                     offset=0.20, width=0.05, colnames = F) +
  theme(legend.position = c(1.05,0.50),
        legend.title = element_text(size = 18), 
        legend.text = element_text(size = 18))+
  scale_fill_viridis_d(name = "Gregariousness (GG)", 
                       direction = -1,
                       breaks = c("solitary", "shoaling","schooling"),
                       limits = c("solitary", "shoaling","schooling"),
                       na.translate = TRUE,
                       option = "D", #'magma',
                       begin = 0.2,end = 0.8,
                       guide = guide_legend(order = 1), 
                       na.value = 'grey')+
  scale_y_continuous(expand = c(0,4))+
  annotate('text', x = 1.23, y = -7.5, label = 'GG', angle = -85, size = 4)

#Check graph
habitat_traits_prob_agg_final

SAVING OUTPUT

There are numerous issues with saving output either via ggsave() or dev.copy2pdf() functions. I’m batching the output saves below and doing all at once.

Basic phylogeny

Including code for saving this graph. Code for saving output for all other figures below are included in the .Rmd and not in the rendered document.

Phylogeny by Class

Phylogeny by diet use

Phylogeny by habitat Use Traits